%% Check coverage of CIs & Compute empirical coverage probability

function [bar_tau, predict_se, empirical_cov, cov_p_value, cov_abs_dev] = check_CI_robust(hat_theta, hat_vartheta, rep_base_se2, rep_val_se2, nu, cov_stat_sim)
% INPUT: estimates hat_theta, hat_vartheta, se^2 rep_base_se2, rep_val_se2,
% EB estimate nu, simulated sample of coverage statistics cov_stat_sim
% OUTPUT: standard error of the predictive distribution predict_se,
% empirical coverage probability empirical_cov, p value of coverage prob

    N = size(hat_theta,1);
    J = size(hat_theta,2);

    % Compute VA and VB for each group i
    inv_terms   = 1 ./ (nu + rep_base_se2);            
    w           = inv_terms ./ sum(inv_terms, 2);  
    w_sq        = w.^2;                             
    VA = (1 + sum(w_sq, 2)) .* nu;                
    VB = rep_val_se2 + sum(w_sq .* rep_base_se2, 2);
    ratio = VA./VB;

    % Compute robust cv
    robust_cv = zeros(N,1);
    for i = 1:N
        mu = ratio(i);
        robust_cv(i) = cva(mu, [], 0.2);    
    end

    % Compute the upper and lower bounds
    bar_tau = sum(w .* hat_theta, 2);

    lower_bound = bar_tau - robust_cv .* sqrt(VB);
    upper_bound = bar_tau + robust_cv .* sqrt(VB);
    predict_se = robust_cv .* sqrt(VB) / 1.282;

    % Check whether the reported estimates lie in the predicted CIs
    hit = (hat_vartheta >= lower_bound) & (hat_vartheta <= upper_bound);

    % Compute empirical coverage prob
    empirical_cov = sum(hit)/N;
    cov_abs_dev = abs(empirical_cov - 0.8);

    % Compute p-value for coverage statistics
    cov_stat_obs = (empirical_cov - 0.8)^2;
    cov_p_value = mean(cov_stat_sim >= cov_stat_obs);

end